﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary
{
    /// <summary>
    /// This class provides the approximation of the jacobian matrix by using extrapolation.
    /// </summary>
    [Serializable]
    public class RealJacobianExtrapolationApproximator
    {
        /// <summary>
        /// Internal array for extrapolation.
        /// </summary>
        private readonly double[] STANDARD_H = {1e-1, 1e-2, 1e-3, 1e-4, 1e-5};

        /// <summary>
        /// Internal array for extrapolation.
        /// </summary>
        private readonly double[] STANDARD_H_MULTIPLY_2 = {2e-1, 2e-2, 2e-3, 2e-4, 2e-5};

        /// <summary>
        /// Internal array for extrapolation.
        /// </summary>
        private readonly double[] STANDARD_H_SQUARE = {1e-2, 1e-4, 1e-6, 1e-8, 1e-10};

        /// <summary>
        /// The source function for the to approximate partial derivative.
        /// </summary>
        private IHardMultivariateRealFunction sourceFunction;

        /// <summary>
        /// Initializes a new instance of the <see cref="RealJacobianExtrapolationApproximator"/> class.
        /// </summary>
        /// <param name="sourceFunction">The source function for the to approximate partial derivative.</param>
        public RealJacobianExtrapolationApproximator(HardMultivariateRealFunction sourceFunction)
        {
            if (sourceFunction == (HardMultivariateRealFunction) null)
            {
                throw new ArgumentNullException("sourceFunction");
            }

            if (sourceFunction.HardRealFunctionPointer == (HardMultivariateRealFunctionPointer) null)
            {
                throw new ArgumentNullException("sourceFunction");
            }

            this.sourceFunction = sourceFunction;
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="RealJacobianExtrapolationApproximator"/> class.
        /// </summary>
        /// <param name="sourceFunction">The source function for the to approximate partial derivative.</param>
        public RealJacobianExtrapolationApproximator(IHardMultivariateRealFunction sourceFunction)
        {
            if (sourceFunction == (HardMultivariateRealFunction) null)
            {
                throw new ArgumentNullException("sourceFunction");
            }

            this.sourceFunction = sourceFunction;
        }

        /// <summary>
        /// Gets or sets the source function for the to approximate partial derivative.
        /// </summary>
        /// <value>The source function for the to approximate partial derivative.</value>
        public IHardMultivariateRealFunction SourceFunction
        {
            get { return sourceFunction; }
            set { sourceFunction = value; }
        }

        /// <summary>
        /// Approximates the specified real function derivatives by using extrapolation and returns the jacobian matrix.
        /// </summary>
        /// <param name="x">The vector x at which the function derivative should be calculate.</param>
        /// <returns>
        /// The specified slope at the positions of x.
        /// </returns>
        public GeneralVector Approximate(GeneralVector x)
        {
            double[] f1h = new double[5];
            GeneralVector tempuri = x.Copy();
            GeneralVector result = new GeneralVector(x.Count);

            for (int i = 0; i < x.Count; i++)
            {
                for (int j = 0; j < 5; j++)
                {
                    tempuri[i] = x[i] + STANDARD_H[j];

                    double vecPlusH = this.sourceFunction.SolveAt(tempuri);

                    tempuri[i] = x[i] - STANDARD_H[j];

                    double vecMinusH = this.sourceFunction.SolveAt(tempuri);

                    tempuri[i] = x[i];
                    f1h[j] = (vecPlusH - vecMinusH)/STANDARD_H_MULTIPLY_2[j];
                }

                result[i] = NevillePolynomialInterpolation.NevilleInterpolation(5, STANDARD_H_SQUARE, f1h);
            }

            return result;
        }
    }
}